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Abstract 

The number partitioning problem consists of partitioning a sequence of 
positive numbers {ai, 02, • • • , a^} into two disjoint sets, A and £>, such that 
the absolute value of the difference of the sums of aj over the two sets is min- 
imized. We use statistical mechanics tools to study analytically the Linear 
Programming relaxation of this NP-complete integer programming. In par- 
ticular, we calculate the probability distribution of the difference between the 
cardinalities of A and B and show that this difference is not self-averaging. 

PACS: 87.10.+e; 64.60.Cn 
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1 Introduction 

Although most of the statistical mechanics analyses of stochastic versions 
of combinatorial optimization problems have focused mainly on the calcu- 
lation of the average cost of the global optima jjj], the tools of equilibrium 
statistical mechanics can also be used to evaluate the average performance 
of simple heuristics as well as that of relaxed versions of the original problem 
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[0. In this paper we study both numerically and analytically the Linear Pro- 
gramming (LP) relaxation of a classical NP-complete integer programming 
problem, namely, the number partitioning problem || 

The number partitioning problem (NPP) is stated as follows Given 
a sequence of positive numbers {ai,d2, . . . ,ajy}, the NPP consists of parti- 
tioning them into two disjoint sets A and B such that the difference 



E a i - E 



(i) 



is minimized. Alternatively, we can search for the Ising spin configurations 
s = (si, . . . , sjv) that minimize the energy or cost function 



E(s) 



N 



(2) 



where Sj 



1 if a., G A and Sj 



-1 if a,j € B. Also of interest is the problem 
of constrained partitions, in which the difference between the cardinalities 
of sets A and B is fixed, i.e., 



N 

E s j 

7=1 



(3) 



so that the cardinality of the largest set is N (1 + m) /2. Henceforth we will 
restrict our analysis to the case where the cu's are statistically independent 
random variables uniformly distributed in the unit interval. 

The interest in the NPP stems mainly from the remarkable failure of 
the stochastic heuristic simulated annealing to find good solutions to it, 
as compared with the solutions found by deterministic heuristics B. The 
reason for that failure is probably due to the existence of order of 2^ local 
minima whose energies are of order of 1 /N Q , which undermines the usual 
strategy of exploring the space of configurations {s} through single spin 
flips. It is interesting to note that a very simple deterministic heuristic, the 
differencing method of Karmarkar and Karp, can find with high probability 
solutions whose energies are of order of l/N alo % N for some a > |J. 
For large N, however, the energies of the solutions found by the differencing 
method are orders of magnitude higher than those predicted by theoretical 
analyses, which indicate that the average global optimal energy Eq is of order 
of y/N 2~ N for unconstrained partitions H, §|. A recent exact calculation 
of this quantity yielded Eq = ^2irN/3 2~ N Q. It must be noted that, in 
contrast with combinatorial problems for which the global optimal energy is 



2 



extensive p]], for the NPP this energy is not self-averaging |j| and hence 
Eq cannot be viewed as a realization independent minimal energy. 

In the LP relaxation we relax the integrality requirement on the Ising 
variables s, so that they become real variables, i.e., Sj G (— 00,00). In order 
to keep these variables finite we impose a spherical constraint on the norm 
of the solutions, 

N 

Y,*i= N - ( 4 ) 

i 

Obviously, minimizing the square of the cost (|2|) with Sj real but constrained 
to obey the condition (Q) yields a lower bound to the optimal (square) cost 
of the corresponding integer programming problem. Moreover, a simple 
gradient descent dynamics suffices to attain those bounds numerically. In 
fact, using a Lagrangian multiplier to handle the constraint (Q) we find that 
the following dynamics minimizes the NPP energy, 

= " i = l,-..,N (5) 

where £ = J2j a j s j an d f] is an arbitrarily small parameter that determines 
the step-size of the descent. 

The remainder of this paper is organized as follows. In section 2 we 
show that the LP relaxation yields a trivial lower bound (i.e. the LP cost is 
zero) for unconstrained partitions. That analysis yields, nonetheless, some 
interesting pieces of information as, for instance, the average energy E c ob- 
tained by clipping (i.e. taking the sign of) the spins of the global optimal 
configurations of the LP relaxation. The average performance of the clipping 
heuristic is studied in section 3, where it is shown that E c tends to the aver- 
age energy of a randomly chosen Ising spin configuration, E c — > y/2N/3n. In 
section 4 we calculate the probability distribution of the difference between 
the cardinalities of A and B for the LP global optimal configurations, V c (m), 
and show that m is not self-averaging. Finally, in section 5 we present some 
concluding remarks. 

2 Linear Programming relaxation 

In the canonical ensemble formalism of the statistical mechanics the average 
value of the optimal energy for unconstrained partitions is given by 

E u = -limT (lnZ u ) a , (6) 
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where Z U (T) is the partition function 



Zu(T) = Y[ [°° dsi S\N-J2 

• J —CO \ 



exp 



E(s) 



(7) 



with E (s) given by Eq. (||). Here <5(x) is the Dirac delta and T is the 
temperature. The notation (. . .) a stands for the average over the random 
variables aj. The limit T — > in Eq. (||) ensures that only the states that 
minimize E (s) will contribute to Z u . We now proceed with the explicit 
evaluation of the partition function (|7|) . Using the integral representation of 
the Dirac delta function we write 



Z U (T) 



d(d( 



oo J -co 
oo 



n 



2tt 
dsj exp 



dx 

-co 2"7T 



,i((+ixN-\(\/T 



(8) 



The integrals over Sj and £ can easily be performed yielding 

IttN-3 rco 

Z u(T) = \ ^r-- / dxe lxN {ixf 



AM 2 

f-OO 



JxN (■ \(l-N)/2 



/CO 
d( e~ 
-co 



\C\/T-ixC 2 /M 2 



(9) 



where M2 = o3- At this stage the integral over £ can be readily carried 
out by assuming T — > 0. The result is simply 



Z U (T) 



1 7T 



N-3 rco 



4M, 



dx (ixf 2 e NG ^ 



where 



G u (x) = ix + — - In (ix) 



(10) 



111) 



In the limit of large N, the integral over x can be evaluated using the saddle- 
point method. Noting that the saddle-point is the imaginary x s = l/2i, the 
unconstrained partition function is finally written as 



Z U (T) 



1 



2tt 2 NM 2 



(1 + ln2vr) 



N/2 



(12) 



Since Z u decreases linearly with decreasing T, Eq. Q yields E u = 0. This 
result was verified numerically using the gradient descent dynamics @ to- 
gether with an adaptive prescription to decrease the step-size ry during the 
descent. 
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3 Clipping heuristic 



An easy-to-implement procedure to generate Ising solutions from the LP 
relaxation solutions is to take the sign of the relaxed spins. The average cost 
associated to this clipping procedure is given by 



E r = lim 



3 



aj sign (sj 



(13) 



where (• ■ -}r stands for a thermal average taken with the Gibbs probability 
distribution, i.e., exp [— E(s)/T] jZ u . The zero-temperature limit ensures 
that only configurations that minimize the relaxed cost @ will contribute 
to this average. To evaluate Eq. ([0]) we introduce the auxiliary energy 



E dip {s) = E(s) + h 
with E(s) given by Eq. (0). Hence 



J2 aj sign ( Sj 



(14) 



E„ 



lim T 



d(\nZ c i ip ) c 



dh 



h=0 



(15) 



where Z c i ip is the partition function ([?]) with E replaced by E c i ip . Introducing 
the auxiliary parameter v = J2j a j sign (sj) through a Dirac delta function, 
the calculation of Z c i ip becomes analogous to that presented before, and so 
we will present the final results only. We find 



E r 



1 

2^ 



/oo poo I ' 

dv\v\ dv e*" ( Y[cos(v aj ) 



/ dv \v\ 

2ir 



dv e l 



SIM) 



(16) 



Assuming that E c does not increase linearly with N, in the thermodynamic 
limit only the regions close to the origin {v = 0) will contribute to the 
integral over v in Eq. (|l~6|). Hence using s\nv/v ~ —v 2 /6 yields 



E r . 




(17) 



We have found a remarkably good agreement between this theoretical pre- 
diction and the properly averaged cost obtained by clipping the spherical 
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spins in LP minima generated by the gradient descent dynamics (|5|). In- 
terestingly, for large N the cost (|l7|) is identical to the average energy of a 
randomly chosen Ising configuration s, defined by 



(18) 



Si = ±l 

which thus demonstrate the complete failure of the clipping heuristic. 



4 Probability distribution of cardinalities 



As the distinct LP global minima will have, in general, different cardinalities, 
in this section we calculate analytically the probability distribution of the 
cardinalities difference defined by 



VAfn) = hm N 
v ; T->0 



Nm 



E s i 



(19) 



Using the definition of the thermal average this equation is rewritten as 



V c (m)=limN(^) a , 

1 — >U Z/„. 



where 



Zm = Y[J°° d Sl 5 (iV-E s l) ex P 



E(s) 



Nm 



(20) 



(21) 



and Z u is given by Eq. (^) . The calculation of Z m is a little more complicated 
than that of Z u , as it involves the evaluation of an additional integral due to 
the extra delta function. However, since the steps are essentially the same 
in both calculations we will present the final result only. In the limits T — > 
and N — > oo we find 



z -L 



1 



N V ir s V 



(l + ln27r) Ar/2 exp 



NM 



2V 



2 2 

m 



(22) 



where V = J2j a "j ~ N J2j a j) • I n the limit N — > oo we use the self- 
averaging property, 

/(°i) = fidaf(a) (23) 
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for any function /, to write M 2 /N = 1/3 and V/N = 1/12 so that Eq. @ 
becomes 

Vc{m) = y^ e - 22V ™ 2 m > 0. (24) 

Hence, the mean is (m) = 1/V2Ntt and the variance, cr^ = (tt — 2) /AttN. 
An important quantity is the ratio r m = y / o^/(m), whose vanishing deter- 
mines the self-averageness of the random variable m. In figures 1(a) and 1(6) 
we present the results of numerical experiments to estimate the dependence 
on N of (m) and r m , respectively, for three types of configurations: (i) the 
global minima of the original NPP obtained through the exhaustive search in 
the Ising configuration space for TV < 26; (ii) the legal, Ising configurations 
obtained with the differencing method (we refer the reader to ref. JjJ for a 
clear presentation of this heuristic); and (Hi) the global minima of the LP re- 
laxation obtained with the dynamics (||). We note the very good agreement 
between the latter estimate and the analytical predictions. In all cases, the 
mean (m) decreases like A r_1//2 as N increases, while r m tends to a nonzero 
value (r m = y/ir/2 — 1 ~ 0.755 for the LP relaxation), indicating that m is 
not self-averaging even in the large limit, i.e., the values of m associated 
to the configurations under study depend on the specific realization of the 
set of random variables {cij}. 

5 Conclusion 

In this paper we have illustrated the usefulness of equilibrium statistical 
mechanics tools to investigate analytically the average performance of stan- 
dard relaxation procedures to generate lower bounds to integer programming 
problems, as well as to characterize specific properties of the minima. The 
failures of the LP relaxation and the clipping heuristic to produce relevant 
results for the NPP yield additional evidence to the extreme difficulty of 
devising heuristics to find near-optimal solutions to that problem. 
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Figure caption 



Fig. 1 (a) Average cardinalities difference as a function of l/N 1 / 2 and (b) 
ratio between the standard deviation and the average cardinalities difference 
as a function of N. The convention is O (LP relaxation), y (exhaustive 
search) and * (differencing method). The solid curves are the theoretical 
predictions for the LP relaxation. 
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